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We review the basics of the coupled-cluster expansion formalism for numerical solutions of the 
many-body problem, and we outline the principles of an approach directed towards an adequate 
inclusion of continuum effects in the associated single-energy spectrum. We illustrate our findings 
by considering the simple case of a single-particle quantum mechanics problem. 
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I. INTRODUCTION 

One of the major endeavors in modern nuclear physics, 
is the ongoing quest for an effective nuclear many-body 
theory, a key element in the attempt to extrapolate ex- 
isting experimental data to regimes which arc not cur- 
rently accessible in the laboratory. Ideally, one would 
like to begin with the experimental data available from 
nucleon-nucleon (AW) scattering experiments J Sl Band 
determine the underlying AW interaction |E IS 111111 II] • 
Subsequent nuclear structure calculations are then nec- 
essary to determine a three- nucleon interaction , to ac- 
count for the remaining differences between theory and 
experiment. 

Direct comparison of theory (calculations) to experi- 
mental data is ambiguous if the nuclear many-body prob- 
lem is not solved accurately. Results of numerical calcu- 
lations can be used to constrain fundamental aspects of 
the theory, provided that one has good approximations 
of the exact result in order to avoid introducing a bias 
due to the approximations involved. Much progress has 
been achieved recently by the Green's Function Monte 
Carlo (GFMC) [T^j collaboration in obtaining accurate 
description of nuclear structure properties of light nuclei 
(A < 10) [Tll[T^| . Such ab initio calculations, based on 
realistic interaction models are key for understanding in 
a quantitative manner the experimental data provided 
by modern accelerator facilities at Jefferson Laboratory 
(JLAB) and the future Rare Isotope Accelerator (RIA). 
We believe that a crucial undertaking in modern physics 
involves the extension of this work to heavy nuclei. 

The limitations of the GFMC method to carry out cal- 
culations for arbitrary size nuclei, originate both from 
the inherent Fermi sign problem associated with a non- 
relativistic system of nucleons obeying Fermi statistics, 
and the fact that in the present implementation of the 
GFMC one carries out explicit summations over the spin- 
isospin degrees of freedom. As such, the relevant num- 
ber of degrees of freedom increases as 2 A (^), and the 
size of the spin-isospin vector space quickly becomes pro- 
hibitive. A possible way around this technical problem, 
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is offered by the Auxiliary Field Diffusion Monte Carlo 
(AFDMC) method where one attempts to sample 
the spin-isospin degrees of freedom in addition to the 
spatial coordinates of the nucleons. This method is cur- 
rently under development, and encouragin g re sults have 
been recently reported for neutron matter |l4| . 

For nuclei with A<4, the nonrelativistic Schrodinger 
equation with realistic nuclear interactions can be solved 
very precisely: In a recent benchmark [l5j carried out for 
a 4 He-like nucleus using the Argonne v 8 ' interaction (cen- 
tral, spin-spin, tensor, spin-orbit and the corresponding 
spin-exchange counterparts), the binding energy calcu- 
lated using seven "exact" many-body formalisms agreed 
within 0.5%. Unfortunately, few many-body formalisms 
are currently capable to offer accurate numerical so- 
lutions of the nuclear many-body problem with realis- 
tic interactions, when heavier systems are concerned. 
For finite nuclei with A > 4, the no core shell model 
(NCSM) [H and the coupled-cluster expansion 

(CCE) are the only methods which are able to directly 
compare results with the GFMC for the same type of 
hamiltonian. Both methods however have serious limita- 
tions related to a slow convergence of the results with the 
size of the model space. As we are about to describe in 
this paper, the shortcomings of the CCE arc not inher- 
ent in nature, and we like to believe that the approach to 
solving the CCE we propose here, will allow the CCE to 
provide accurate theoretical output for a medium to large 
domain of masses, in the framework of realistic nuclear 
interaction 0, 0, [hJ and nuclear currents [2(| • We call 
this new approach the continuum CCE (c-CCE), because 
continuum effects are accurately taken in account. 

The coupled-cluster expansion (CCE), also called the 
cxp(S) method, was developed in the early 1960s by Co- 
ester and Kiimmel 0, . While the method is viewed 
as exact, approximations are introduced stemming from 
truncations in the CCE equations, as well as truncations 
in the model space. Practical approaches for nuclear 
structure applications have been notoriously difficult to 
realize. It was not until the 1970s that Zabolitsky and 
Kiimmel [2j| were able to carry out the first detailed cal- 
culations for finite nuclei, using a representation of the 
wave function in coordinate space together with common 
interactions of the time. The results for the binding en- 
ergy of nuclei such as 4 He, 16 O and 40 Ca, were well above 
the Cocstcr line and were taken as evidence for the pres- 
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ence of 3iV and higher-order interactions. 

While the CCE was extensively used in other areas in 
physics and chemistry (see Refs. (2J,[25| for additional in- 
formation), further applications to nuclear structure cal- 
culations based on realistic nuclear interactions proved 
difficult to achieve. In the 1990s we mainly note the at- 
tempt of constructin g a translationally invariant CCE in 
coordinate space [26l . l27l . l28l ]. method which might yet of- 
fer an attractive alternative for finite nuclei calculations 
(see Refs. [2|| [3(j for information on recent progress on 
this approach.) To date, however, these attempts are 
confined to semi-realistic interactions with a V4 structure 
(scalar plus spin-, isospin-, and space-exchange terms), 
and the wave functions are not sophisticated enough to 
provide the realistic account of nuclear forces required in 
the context of realistic nuclear interactions, such as the 
Argonne vis and corresponding thrcc-nuclcon forces. 

Motivated by the availability of more sophisticated 
NN interactions, and riding the wave of the ongo- 
ing expansion in computer power, Heisenberg and Mi- 
liaria |3ll I32I I33I ] reexamined the coupled-cluster ap- 
proach and applied it to the spherical nucleus 16 0, to- 
gether with a realistic nuclear interaction and currents 
which were purposely consistent with the work done by 
the GFMC collaboration: Calculations based on a real- 
istic two-body hamiltonian (Argonne V14 and vis) were 
reported in Ref. [3l|], while effects of a phenomenolog- 
ical three-body interaction (Urbana IX) were reported 
in Ref. [321 ]. An extension of the later application, in- 
volved the microscopic calculation of the charge form 
factor in 16 O. Results were presented in Ref. [33j, and 
showed good agreement with the available experimen- 
tal data [34| when both three-body forces and meson- 
exchange currents were taken into account. 

In its present implementation the CCE equations are 
solved for the special case of spin-isospin shell-saturated 
nuclei, such as 4 He, 12,14 C, 14,16 0. In order to depart 
from the closed-shell case one resorts to an equation of 
motion approach, using the ground-state of the closed- 
shell nucleus as a new vacuum. As such one can describe 
the properties of hole-state nuclei such as 13 C and 15 N, 
and preliminary results are promising. However, such 
an approach will inevitably compound the effect of new 
and old approximations. Also, the present approach has 
limitations due to the intrinsic cutoffs one has to impose 
when defining the model configuration space. The single 
particle states are expanded in a harmonic oscillator ba- 
sis, and satisfy the same type of boundary conditions for 
both hole and particle states. In effect we discretize the 
continuum part of the self-consistent one-body mean-field 
hamiltonian used to define the single-energy spectrum. 
This results in the necessity of a large 50fi.fi configura- 
tion space and subsequent significant storage problems 
and lengthy execution time. 

One would like to obtain exact solutions for the 2-, 3- 
and 4-body systems via CCE. Since contributions due 
to n > A particle-hole configurations are always can- 
celled exactly, the above few-body systems will provide 



scenarios when the truncation of the CCE hierarchy for 
A = 2, 3 and 4, respectively, becomes exact, and we will 
be able to reliably ascertain the numerical accuracy by 
direct comparison with similar results obtained via other 
exact methods. Moreover, by solving the two- and three- 
body problem, we will effectively eliminate our present 
reliance on the closed-shell hypothesis, and allow for di- 
rect calculations of nuclei inside the shell via CCE. Given 
our past experience regarding the ground state calcula- 
tion of spin-isospin shell-saturated nuclei in configuration 
space, we can expect that a four-cluster truncation of the 
CCE equations will allow an accurate description of ar- 
bitrary nuclei beyond A = 16. 

From a numerical standpoint, we must always remem- 
ber that we have to solve a finite physical system, using 
finite computational resources, in a finite amount of time. 
We submit that the resolution of the many-body problem 
is linked to the successful design of a numerical algorithm 
that allows for an efficient implementation on a massively 
multiprocessor machines. In order to preserve the scala- 
bility of a numerical algorithm on a multiprocessor ma- 
chine, we propose to address the problem of finding the 
energy spectrum of the hamiltonian by strictly following 
the two commandments: 

1. Thou shalt not solve matrix eigenvalue equations. 

2. Thou shalt not perform explicit matrix inversions. 

In order to solve this conundrum, we must pursue those 
iterative numerical schemes which are most likely ap- 
proachable using Monte Carlo techniques. It is our con- 
tention that the basic expertise for solving the nuclear 
many-body problem is available at this time: The GFMC 
is the only many-body method available today that obeys 
the above constraints. However, in its present formula- 
tion, the GFMC carries out explicit summations over all 
possible spin-isospin combinations, and this leads to a 
2 A (^) increase for the number of the relevant degrees 
of freedom. We propose to investigate the possibility of 
solving a comparatively small (n < 4)-cluster truncation 
of the CCE hierarchy of equations using Monte Carlo 
techniques. Assuming the CCE truly exhibits a 1 /A- type 
convergence, then it is reasonable to expect that approx- 
imations based on the CCE will accurately describe nu- 
clear structure with realistic nuclear interactions and cur- 
rents, around and beyond 16 0. The CCE equations will 
still involve explicit spin-isospin summations, but these 
will entail only minimal complications because the small 
size of the cluster will drastically simplify the problem. 

To summarize, our long-range interest is two-fold: 
firstly, we want to develop a continuum version of the 
CCE in an attempt to efficiently describe physics at large 
distances; secondly, we are interested in finding an effi- 
cient numerical procedure to solve the CCE equations. 
We find it useful to begin our investigations with the 
simplest dimensional realization of the many-body prob- 
lem: the one-body case, A = 1. The mere simplicity of 
the problem will allow for an exhaustive, detailed study 
of the CCE equations. 
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II. SINGLE-ENERGY SPECTRUM 

Consider the Hilbert space of states associated with 
the one-body hamiltonian 



H = T + V 



(1) 



and a basis set in this Hilbert space. In general, this basis 
features both a discrete and a continuum part, such that 
any state in the Hilbert space can be expanded as 



d 3 p 



(2) 



That is to say that the basis states in the spectrum satisfy 
the orthonormality relations 

(4>n\(t>m) = Sum , (3a) 

(fa\<j> p ) = (2tt) 3 5(p—p') = S pp , (3b) 
(&,|&r) =0, (3c) 



and the basis set is complete 



^2 IMi^n 



d 3 p 



|#)(#|= 6(r-r'). (4) 



The above basis can be the same as, but is not re- 
stricted to, the set of eigenvectors corresponding to the 
single-energy spectrum defined by the solutions of the 
Schrodingcr equation 



Ho 



(5) 



We propose to start with a known set of discrete single 
particle wave functions {4> n } in the above Hilbert space, 
which satisfy Eq. (|3a|l . and subsequently construct a set 
{4>p}, which obeys the remaining orthonormality condi- 
tions, Eqs. (|3bl3c|l . the closure relation Q, and satisfies 
a predetermined boundary condition (i.e. bound state or 
scattering state). 

This can be done as follows: In order to satisfy the 
orthogonality conditions between the appropriate con- 
tinuum single-particle wave functions, we start with a 
complete set of continuum wave functions 



\Xp) = \Xp) 



^2 \M(<f>-n 



\Xp) 



(6) 



For instance, if the continuum wave function Xp(^) is just 

a plane wave, Xp( 7? ) = e lp r , then (4> n \Xp) = 4>n{p) denote 
the Fourier transform of the wave function <j) n (r) . The set 
of continuum wave functions {xp} are orthogonal to the 
discrete spectrum of the hamiltonian, i.e. (4> n \Xp) = 0, 
and satisfy the completeness relation 



^2 \M(^>n 



d 3 p 
(2^)3 



\Xf)(Xp\ = S(r-r') 



(7) 



These wave functions however are not orthogonal to each 
other 



ixp\xp) 



3pp' 



(8) 



Without loss of generality we can orthonormalize this set 
{Xp} to obtain the desired set of continuum wave func- 
tions, {4>p-}- Note that there is no inconsistency here, 
since Xp(f) is not simply a linear superposition of plane 
waves, but we can define </>p(r) as a linear combination 
of Xp{^) m order to perform the orthonormalization. In 
practice, the fact that the definition of Xpif) involves the 
Fourier transform of <^ n (r), establishes the appropriate 
grid representation of the single-particle wave functions. 
The set {<f>p} is obtained by orthonormalization of the set 
{Xp} a t these grid points. This concludes our construc- 
tion. 

The fact we actually have to normalize the continuum 
basis states on a grid, may be a serious drawback. After 
all we started this process with the idea in mind that a 
formal discretization of the continuum prevents us from 
obtaining rigorous numerical convergence with modest 
computational resources. It is therefore unfortunate that 
an actual discretization appears to be necessary at this 
point. However, we claim that already at this point we 
are better off in this formalism since all basis states, and 
most importantly the continuum basis states, have the 
correct boundary conditions built in. Moreover, we will 
see that the above construction is not binding us into a 
certain course of action. We will argue later, that we 
can in fact circumvent these apparent difficulties, and 
return to a representation formulated entirely in coordi- 
nate space where the explicit knowledge of the continuum 
states is not mandatory. 



III. SECOND QUANTIZATION 

We proceed now by outlining the second quantization 
language we will subsequently use. We start by intro- 
ducing the field operators ^(r) and ^(r), as linear com- 
binations of creation and destruction operators c^ct,, 
respectively: 



(9) 



where the single-particle wave functions 4>(r) = 
{</>„(r), (j)j;(r)}, form a complete basis, and satisfy the 
orthonormality properties, as explained in the previous 
section. The inverse transformations are 



bt = J d 3 r C(?W) 
up = / d 3 r 0;-(r) i>{r) . 



(10) 
(11) 
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In other words, removes a nucleon from the hole orbit 
^n(r); while ELp removes a nucleon from a particle orbit 
4>p{f). With these definitions, we can write 

c t = (4< ife P> £ ^> (12) 
[b n , if e n < e F , 

and 

c = luf il £p - £F > (13) 
[bT , if e n < e F . 

Here, Ef denotes the Fermi energy. Note that in the re- 
gion between the Fermi energy and the separation energy 
(e p = 0), the "continuum" spectrum, |</y), is actually dis- 
crete, and the integral has to be interpreted as a sum. 

Since we are interested in many-body systems obeying 
fermi statistics, we ask the field operators to obey the 
anti-commutation relations 

{^(r), $(r*)} = 5(r - ?) . (14) 

Spin and isospin indices are implied. In turn, the creation 
and annihilation operators obey the canonical relations: 

{ C L c m\ = S„m , {ct, Cp } = , 

{c„, c m } = {4, c+J = {c p -, Cj? } = {ct, ct} = , 
{c„, c p -} = {c^, el} = {cl, Cp-} = {c„, ct} = , 

or 

{bjj,b m } = S nm , {at, a^/} = 8~p , 
{b„, b m } = {bl, bjj = {ap-, a p -,} = {at, at,} = , 
{b„, a p -} = {bt , at} = {b* , a p -} = {b„, at} = . 

More definitions are in order: We introduce the bare 
vacuum, |0), such that 

$(r) |0) = , (17) 

or 

ap-|0)=0, bt|0)=0. (18) 

In turn, the physical vacuum, |$), which will play the 
role of the reference state for exp(S), is defined as the 
exhaustive collection of minimal configurations of a given 
symmetry, obtainable from the bare vacuum: 

a^|$) = 0, b„|$)=0. (19) 

We have 

41$) = 4|$) , c„|$) = bjj$> . (20) 

Finally, the second quantization representation of the 
one-body hamiltonian operator is given by 

H = ^IHI^c/,. (21) 

a/3 



By the same token, we will represent for instance, the 
Iplh cluster-correlation operator, as 

Si = E /(04%" b ^- ( 22 ) 



IV. CCE EQUATIONS 

In a traditional shell-model approach, one calculates 
the matrix elements of the hamiltonian 

(*x|H|* 2 ) = E^^fanlHl^,) (23) 

mn 

+ E/ (§3 ^(^IHI^) 




Then, for an arbitrary eigenstate of the hamiltonian, 
we write the Schrodingcr equation 

H E l*a)(*a|*) = ^ E . ( 2 4) 

a a 

which leads to the usual eigenvalue problem 

E(<*«|H|*/j)-£?„<ya/j)<*/j|*) - 0. (25) 

We make the convention that Greek indices run over both 
the discrete and continuum spectra, i.e. the sum J2 a ^ s 
in fact a sum over the discrete, and an integral over the 
continuum. When solving the above eigenvalue, one is 
usually forced to discretize the continuum: one uses a 
discrete basis, which results into a finite matrix, which 
can then be diagonalized using regular linear algebra nu- 
merical packages. 

We would like to avoid this approach. Therefore, 
rather than calculating the "entire" spectrum of the 
hamiltonian at the same time, we will attempt to cal- 
culate the spectrum of the hamiltonian one state at a 
time, similar to the GFMC approach. In the same spirit, 
we will attempt to calculate observables without explicit 
knowledge of the actual wave functions: wave functions 
are not observables, and we will concentrate on develop- 
ing a formalism allowing direct calculation of the expec- 
tation values of operators. 

In the CCE formalisms, one first introduces an ansatz 
for the eigenstate \$f) of the hamiltonian 

I*) = e s |$) , (26) 

where |$) is the physical vacuum or reference state. The 
ansatz for \^f) satisfies the normalization condition 

(*|$) = 1. (27) 
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Here S denotes the many-body cluster-correlation oper- 
ator: 



A 

E 

n=l 



(28) 



where denotes an npn/i-configuration creation oper- 
ator. Mathematically speaking, the interpretation of a 
npnh configuration is that an element of a basis set span- 
ning the many-body Hilbert space. For a given npnh con- 
figuration, the correlation function for n-nucleons, or the 
npnh correlation coefficients S n , is simply the associated 
expansion coefficient. 

We emphasize that the physical vacuum, |<I>), is de- 
fined as the minimal set of configurations obtained from 
the bare vacuum, which obey a certain set of symmetries 
and boundary conditions. Appropriate choices can be 
made for |$) which will allow us to independently con- 
struct both the ground state and excited state spectrum 
of the hamiltonian, without changing the single-energy 
spectrum. The cluster-correlation operator is a rank-zero 
tensor operator and will leave unchanged both the sym- 
metry and boundary conditions of |$). Therefore, the 
information regarding a particular set of quantum num- 
bers is contained in the physical vacuum |<f>); properties 
such as orthogonality of the various eigenstates will be 
satisfied by construction. In so doing, we depart from 
prior attempts to obtain the excited state spectrum of 
the hamiltonian by means of an equation of motion ap- 
proach which relies on a prior calculation of the ground 
state. 

With these definitions, the Schrodinger equation for 
the eigenstate \^/) of the hamiltonian corresponding to a 
chosen symmetry, 



hi*) = E\m 



becomes 



or, in normal ordered form, 
(e- s He 



l*> = E |<f>> 



(29) 



(30) 



(31) 



where the subscript c indicates the pure creation part of 
a normal ordered operator. In principle, the resulting 
system of integral equations must be solved for the npnh 
correlation amplitudes S n and the energy E. It is desir- 
able, and in the spirit of the approach we initiate in this 
paper, to design a way of calculating E without explicitly 
calculating (and storing!) S n . 



V. OBSERVABLES 

Consider the expectation value of an arbitrary opera- 
tor O 



(32) 



We perform the following transformations 
_ ($|e st Oe s |$) 

Jm¥) 

($|e st e s e- s 0e s |^) 

We can insert a complete set of npnh configurations, and 
obtain 

= £ W<*i°.«-wi»>. 

n VI/ 

where in the last step we have used the fact that 
and S commute. We introduce a new cluster-correlation 
operator, S, by its many-body decomposition, as 

'-E^-EU. (33. 



Note that 



(*|Ot|tt) (*|O n |*) 



(34) 



and therefore the S n amplitudes are real. With this def- 
inition, we can calculate the expectation value of the op- 
erator O as 

o = ($|S e- s Oe s |$) = ($|e st 0e- st S f |$) . (35) 
Based on the (second) definition of the npnh amplitudes 



Sn 



($|S e" s O, 



(36) 



we can determine S n via an iterative procedure. 

Note that the correlated ground state |*) is not a 
translational-invariant wave function. Therefore, in prac- 
tice one must correct for the effects of the center-of-mass 
(CM) motion. In our previous studies, a many-body 
expansion has been devised to evaluate the corrections 
required by the calculation of observables in the CM 
frame [35l | . The accuracy of the proposed procedure was 
tied to the success of a good separation of the CM and 
relative coordinates degrees of freedom in the CCE solu- 
tion. This hypothesis was tested [32] by adding to the 
Hamiltonian a purely CM piece, multiplied by a strength 
parameter, and making sure that the binding energy is 
independent of this parameter. 

We review this issue now on fundamental grounds: In 
order to produce translational-invariant many-body wave 
functions, one must seek a simultaneous eigenvalue of 
both the hamiltonian and the total momentum of the 
system: 



H |*) = E |*) 

p m = P i*) 



(37) 
(38) 
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or 

e" s He s |$) = E\§), (39) 
e -s Pe s = p |$) . (40) 

A closer inspection of the above set of equations shows 
that this not an over-determined system of equations: 
By definition, the hamiltonian is the internal hamilto- 
nian of the many-body system, where the CM kinetic 
energy has been removed. The following commutators 
are identically zero: 

[H, P] = [e" s He s , e- s Pe s ] = 0. (41) 

One can replace Eq. I|4(J|) with a linear combination of 
Eqs. and : 

e" s [H, P]e s |$) =0. (42) 

Since the hamiltonian and the total momentum com- 
mute, it follows that this equation is always satisfied by a 
solution of Eq. (|39[) . Therefore, the dressed hamiltonian, 
e~ s He s , docs not modify the CM degrees of freedom 
content of the physical vacuum, |<J>). 

In practice one has to truncate the CCE hierarchy of 
equations, and an inadequate truncation may render the 
above discussion inconsequential. However, the above 
statements remain valid for a truncated version of the 
CCE, provided that a truncated cluster-correlation oper- 
ator S( w ) can be defined, and satisfies the identities 



and 

e~ sm Oe sW =0+[0,sW] (44) 
+ I[[0 ! SW] ! SW]+---. 

An acceptable truncation of the CCE hierarchy is gener- 
ated by the prescription 

S„ = 0, n > N . (45) 

The physical vacuum plays indeed a key role in the CCE 
formalism: When done correctly, and even if |$) breaks 
certain symmetries of the desired solution, the CCE will 
leave the center-of-mass part of the physical vacuum un- 
changed. Therefore, one can perform corrections at the 
level of the physical vacuum, and the further inclusion 
of the npnh correlations will leave these corrections un- 
changed. This is an important observation, since we have 
the physical vacuum in closed form, and we would rather 
not explicitly calculate the npnh correlations. If |<J>) can 
be written as the product of a CM wave function and a 
wave function of the system with respect to the center- 
of-mass, then one can carry out the CM corrections via a 
simple correction factor, just like in the shell model with 
harmonic oscillator wave functions. 



VI. A = 1 : SINGLE-PARTICLE QUANTUM 
MECHANICS 

So far, our discussion has been quite general. We will 
particularize now to the simplest possible case, namely 
the case of a single-particle quantum mechanics problem. 
The CCE is a many-body formalism, but the procedure 
is quite general, and we can apply it to the A = 1 case. In 
this context, we will continue using the many-body lan- 
guage, and refer to "configurations" and "correlations" , 
but solely because this is an intrinsic part of the CCE 
formalism. 

In the one-body case, the cluster-correlation operator 
is particularly simple 

S = Sj . (46) 

Since the cluster-correlation operator contains only Iplh 
correlations Eq. I|26|) can be written as 

|*) = e s |$) = (1 + S) |$), (47) 

where no matter what symmetry we are interested in, 
there is only one possible configuration in |$). We in- 
troduce the physical vacuum for the one-body case, as 

l*> = \M- (48) 

For simplicity, we will consider that \<ph) is the only state 
in the discrete spectrum satisfying the required symme- 
try. An extension to the case when we have more than 
one discrete state in the spectrum is straightforward and 
the relevant equations will be relegated to Appendix [BJ 
In order to derive the corresponding CCE equations 
we use the identity (see also Appendix^} 

e- s He s = H + [H, Si] + | [[H, Si], Si] . (49) 

Then, Eq. (|31fl is equivalent to the following set of equa- 
tions 

E = (H) + ([H,Si]) , (50) 
= (H) 1 + ([H,S 1 ]) 1 + 1([[H,S 1 ],S 1 ]) 1 . (51) 

These equations are solved for the energy E correspond- 
ing to |\&), and the correlations S\. 

Using Eqs. (|46[) and (|47|) we can write the solution of 
the Schrodingcr equation as 

I*} = \4>n) + J 7^ SffHfa) . (52) 

We note that for the one-body problem, it is particularly 
simple to find an interpretation or an independent check 
for the CCE equations. Based on Eqs. l|T7|) and Q29[l. we 
can readily derive 

E = (^ h |H|*) , (53) 
ESgh = (^ P |H|tt). (54) 



7 



Using the commutators listed in AppcndixlAl the CCE 
equations we need to solve become 



E = 



|H|0 



d 3 p 



b h \H\4>p) Sp h 



= (</> p -|H|^ 
d 3 Pi 

(27T) 

d 3 pi 



(2tt) 3 

3 wvW'^Vi) Spin — Sph {<j>h\K\(t>h) 



(55) 



(2tt; 



3 Spft (^/ilHl^p-j) Sp-^ft 



(56) 



By inspection, the last equation gives 



ESp h =(<p^H\4> h ) 



d 3 Pi 
(2tt) 3 



MH|^)5 plh . (57) 



Equations (|55|) and 157() involve explicit integrals over 
correlations in momentum space, and we will refer to 
this set of equations as the P-set. These equations are 
entirely consistent with Eqs. i|53|) and (|54|l above. 

We note for completeness, that by formally introducing 
the notations 



T 



bh\H.\<h) 



(58) 
(59) 
(60) 



the solution of P-set of equations can be formally ob- 
tained as 

Si = nrv (6i) 

E = (0h|H|0 h ) + v T ■ [A]~ l v. (62) 

This implies, of course, that in order to get the energy, 
one has again to discrctize the continuum, and the solu- 
tion is obtained by inverting a matrix. The procedure is 
iterated until convergence is achieved. We do not propose 
using this approach. 

The above momentum space integrals can be avoided 
in favor of a coordinate representation: We begin by in- 
troducing the /io/e-function 



Xh(r) 



\Xh) 



4>p(r) Sph 



d 3 p 
(2^)3 

(2^ I*?) S ^ 



(63) 



(64) 



Note that in this simple one-body case, we have 

\ Xh ) = |*) - \4>h) ■ (65) 

With this notation, we obtain the equivalent set of equa- 
tions 

£?=(0 fc |H|0/i) + (0h|H|xfc), (66) 
ESp h =(0 # |H|^) + {4>p\H\ Xh ) , (67) 



where the Sph correlations are interpreted as projections 
of the ZioZe-function, \xh)i onto the continuum set of 
states, i.e. 



S,7l 



ph 



{<t>p\Xh) 



(68) 



We can transform Eq. I|67|) into an equation for the hole- 
function: we multiply by \(f>p) , and integrate over the mo- 
mentum variable. Using the closure relationship, Eq. Q), 
we obtain 

E\ Xh ) = H\<j> h )+H\ Xh ) (69) 
-|^){<0 h |H|^) + <^|H|x h )}. 

This equation is consistent with the eigenvalue equation 
for |*) : 

(H - E) \ Xh ) + (H - E) \4> h ) = (H-P)|*) = 0, 

and can be used to evaluate \xh)- We have : 

(H-E)\ Xh ) = -(H-E)\j> h ), (70) 

subject to the constraint that \xh) is orthogonal to the 
discrete basis states, i.e. 



K\xh) = o. 



(71) 



At this point we consider a new complete set of functions 
in coordinate space, which satisfy the required boundary 
conditions. Since these functions may overlap with the 
discrete states of the original basis, we will remove this 
overlap by using the procedure outlined in Eq. ©. We 
will then expand \xh) m the new basis, and use Eq. I|70|) 
to determine the corresponding expansion coefficients. A 
similar technique is successfully used in continuum RPA 
theory [3(|. The approach is also similar to a method 
introduced many years ago by Podolsky [23, [2~3 , which 
allows for the substitution of an integral over the contin- 
uum part of the spectrum in terms of a bound-state eigen- 
function of an auxiliary hamiltonian. Podolsky's method 
was initially introduced to treat dispersion in hydrogen 
atoms, and it has been recently applied to the calcula- 
tion of the electric polarizability of the deuteron (3^ with 
various realistic NN potentials. 

We will refer to the set of equations (|66|l and J7UJ), as 
the P-set. The key element in this representation is the 
fact that an explicit knowledge of the continuum part of 
the basis states, \4>p), is not mandatory. 

We turn our attention now to the calculation of ob- 
servables: For a one-body problem, S contains only Iplh 
configurations 

S = 1 + Si. (72) 
We determine Si as a solution of the equation 



Sp, = (<S>\e- s b h ape s \<S>) 
d 3 pi 



(73) 



r Spih (^Ibha^ e s b ft ,a J? e s |$) 

(ZTT J j 
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and subsequently, observables will be calculated as 



VII. ITERATIVE SOLUTIONS 



o = ($|e- s 0e s |$ 
d 3 p 



(74) 



(2nf 



After working out the necessary commutators, we obtain 



ph 



ph 



(27r) 3 Pl p 



d 3 p 



(75) 
(76) 



These equations involve integrals in momentum space, 
and therefore, Eqs. I|75[l and (|76|l represent a natural 
complement to the P-set of equations discussed above. 
In order to derive the corresponding coordinate repre- 
sentation, we need to introduce a modified ZioZe-function 



Xh(r) 



\Xh) 



d 3 p 
(2^F 



d 3 p 



<t>p(r) Sph 



'p') Sph 



(77) 



(78) 



Note that by construction we have (Xh\<Pn) = 0. With 
this notation, the above Eqs. i|75|) and (|76|l become 



(Xk\ = (Xk\ l--^<Xh|H|*j 



o = (^ h |0|*) + {x h |0|*)--<Xh|H|*)<^|0|* 



(79) 
>■ 

(80) 



These equations allow for the calculation of the observ- 
able o entirely in coordinate space. 

An interesting expression for the expectation value of 
an arbitrary operator is obtained by combining Eqs. I|79|) 
and <(8"0"j) . We obtain 



o = (#101*) 



i-i(x fc |H|*) 



(81) 



This last equation is a direct consequence of the fact that 
CCE does not automatically generate a normalized wave 
function |#). Of course, for the simple case we consider 
here, a brute force normalization of |#) is always possible, 
but such a procedure will be impossible for A > 1. Our 
prescription (|7TJ|> for calculating the expectation value of 
O provides an implicit normalization for the purpose of 
calculating the observable o, by means of a simple sub- 
straction from the expectation value calculated using the 
unnormalized wave function |W), as shown in Eq. I|81|l . 



Let us consider the P-set of equations, Eqs. (|55jl and 
(I57|) . An iterative solution for the lplh correlations, Sph, 
can be obtained by iterating Eq. (|57fl . We have 

S pih = -g; (0ft |H|^) + i (cj) ffl \U\cl)f 2 )((/)p 2 \H.\^ h ) 



1 



* IHI 



» P -JH|^) 



(82) 



Implicit integration over all repeated continuum indices 
is assumed. By combining Eqs. (|82[) and l|55|) . we obtain 
a polynomial equation for the energy of the desired state 



1 



E 



— ^ \n\4>p 2 ) ((j)p 2 \n\4> h ) 



E 



>ft|H|0 fc ) 



yJH|0 h ) + 



(83) 



The correlations Sph do not explicitly enter the last equa- 
tion. 

The above equations are not suitable for a Monte Carlo 
approach to obtaining the energy E, since the energy 
itself appears as a inverse power factor in both Eqs. I|82|) 
and (|83fl . This situation is an artifact of the one-body 
problem we consider here, and can be traced back to 
the simple substitution made in deriving Eq. 1|57|1 from 
(|56[) . For the actual many-body (A > 1) problem, the S n 
equations will also feature an energy denominator, but E 
will be replaced by the npn/i-energy, e n = {ep 1 — en, ■ ■ ■}■ 
The existence of an energy gap between the discrete and 
continuum spectrum will provide a small parameter for 
this expansion. 

Nevertheless, Eq. I|83fl deserves a second look: We use 
the closure relation (0} to eliminate the integrals over the 
continuum basis states. We obtain 



E 



W\<t>h) 



|{(0 h |H a |^>-[(^|H|^>] 2 } 

-2(0 /i |H 2 |0, l )(^|H|^) 
[(0,|H|^)] 3 } + ... (84) 



This equation can be used in practice to calculate E for 
A = 1, by considering successive truncations of the above 
series, and finding the root of the resulting polynomial 
equation in E. 

An iterative scheme is also possible for calculating ob- 
servables. The observable o can be obtained by em- 
ploying Eq. I|79|) and a repeated substitution of (xh\ m 
Eq. dSH>. We have 



00 c— iV" 



k=0 



(85) 
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To conclude, let us describe another traditional strat- 
egy for obtaining the numerical solution of Eq. I)57p. 
which also relies on an iterative procedure: We begin 
by making an initial guess for the \plh correlations: 

^ = -<«M, (86) 

where we use the notation e a = (^ a \H\^) a ) , the diagonal 
matrix elements of the hamiltonian. Then, we improve 
on this initial guess, by finding a first-order correction to 
Sph m t ne standard fashion: 

— Sph + eph ■ (87) 

This procedure is iterated until convergence is achieved. 
The correction eph is the solution of the equation 

E° ep h - J |^ (#|H|^ ) ep lh (88) 

+ Sph J (0ft I H|^) ep lh 

= (^|H|^) + J |^ (0 P HH|^) S* fc ~ E° S% , 

where E° = E[Sp h ]. Of course, this procedure also re- 
lies on discretizing the continuum, and involves a ma- 
trix inversion required for inferring the correction eph as 
solution of the above system of linear equations. How- 
ever, given a particular hamiltonian H, it is important to 
ask the question wether this procedure is stable and con- 
verges to the correct solution. If this is indeed the case, 
then one can show that the solution of the CCE equa- 
tions can be obtained as the sum of an infinite Neumann- 
like series, and one can envision attempting a numerical 
summation of this series by means of a nonlinear Monte 
Carlo approach. The detailed discussion of this approach 
is beyond the scope of the present paper. 



VIII. NUMERICAL DIGRESSIONS 

At this point, a numerical study of some of the above 
formal statements might be helpful to illustrate our point 
of view, and the case of a simple single-particle quantum 
mechanics problem provides interesting insight in the na- 
ture of the approximations. Let us consider a spherically 
symmetric Wood-Saxon (WS) potential well 

V(r) = -V /{l + exp[{r-R)/a]} , (89) 

and address the problem of finding the bound-state spec- 
trum of the associated one-body hamiltonian. The choice 
of a WS potential is not entirely accidental: the WS po- 
tential may play a special role in our future many-body 
studies, as it mimics well the mean-field nuclear poten- 
tial, and accurate numerical methods for calculating the 



bound state energies and wave functions are readily avail- 
able in the literature. As such, we may use single-particle 
WS wave functions to construct good variational repre- 
sentations of the physical vacuum in the CCE. 

We consider a WS potential well with depth Vq = 100 
MeV, radius R = 3 fm, and diffuseness parameter a = 0.6 
fm. In the following we focuss on the calculation of the 
bound-state energies, and the corresponding probability 
densities. Numerical results are collected in Tableland 
Figs. Eh-/. 

We begin by employing two textbook 01 methods for 
finding the bound-state eigenvalues of the hamiltonian. 
Firstly, in order to establish a baseline for our compar- 
ison, we employ a shooting technique and solve the ra- 
dial Schrodinger equation on a grid using the Numerov 
method. The radial wave function calculation begins at 
the origin and the numerical solution is stepped out for 
a fixed energy E, to r max = 12 fm in 2 12 x 10 3 steps. 
In the presence of an eigenvalue, E n i, the wave function 
flips between large negative and positive values for small 
variations in the energy E. The energy E is fine tuned 
in order to capture this transition. In this limit E tends 
to E n z, the true eigenvalue of the Schrodinger equation. 
For E = E n £, the numerical wave function becomes nor- 
malizable over the radial interval [0,r max ], and we can 
compute the associated density probability on the inter- 
val. These eigenvalues and densities will be referred to 
as "exact" for the purposes of our numerical comparison. 

Secondly, we obtain the solution of the Schrodinger 
equation in a harmonic oscillator representation, by solv- 
ing the eigenvalue problem illustrated in Eq. (|25|) . In this 
(Heisenbcrg) picture, the continuum part of the hamil- 
tonian spectrum is inherently discrctized. The harmonic 
oscillator (HO) wave functions depend on the choice of 
the oscillator parameter 6, and this in turn introduces 
a dependence of the calculated spectrum on the oscilla- 
tor parameter b. The corresponding cigenfunctions are 
linear superpositions of HO functions, which partially 
compensates for the fact that the HO functions do not 
exhibit the corrected asymptotic behavior demanded by 
the Wood-Saxon potential. We consider a number of 
N g = 68 functions in our HO expansion, to correspond 
to the N g = 68 collocation points used in our Gauss- 
Hermite quadrature scheme. Test runs for N g = 52 show 
that both the bound-state eigenvalues and eigenvectors 
of this hamiltonian have largely saturated with N g , and 
the subsequent improvements arc only minute with in- 
creasing N g > 52. The eigenvalues are very close to the 
expected exact results. The calculated densities for the 
lower energy states are reasonably close to the exact re- 
sult: the more bound the state is, the better the density 
agreement. For the shallow (E > — 10 MeV) bound- 
states however, the calculated densities compare poorly 
with the exact results: the model space is simply too 
small to allow for a good representation of the wave func- 
tion, even though the eigenvalue calculation has arguably 
converged already for this space size. 

We turn our attention now to the CCE solution. In Ap- 
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TABLE I: Energies [MeV] of WS bound states: calculated values using the shooting method (exact), direct matrix diagonal- 
ization in a TLO basis (Heisenberg), physical vacuum energy (Eh), and CCE using a pure TLO basis (exp(S): TLO), and a 
combination of a TLO and £ basis (exp(S): TLO), respectively. Where TLO basis are concerned, results are reported for three 
values of the oscillator parameter 6, centered around the optimum b value for the variational ground state energy, b — 1.4 fm . 


state 


Heisenberg 

b = 1.3 fm 


E h 


exp(S): TLO 


exp(S): TLO + C 


exact 


[£ = 0,n = l] 


-71.604115 


-71.154641 


-71.604115 


-71.604093 




[£ = 0,n = 2] 


-27.101582 


-26.341419 


-27.487394 


-27.101544 




[£ = 1, n = 1] 


-50.549663 


-50.189162 


-50.549663 


-50.549664 




[1 = 1, n = 2] 


-7.974057 


-4.649952 


-8.272408 


-7.974058 




[£ = 2,n = l] 


-28.555288 


-28.070778 


-28.555288 


-28.555288 




[£ = 3, n = 1] 


-6.890894 

b = 1.4 fm 


-5.580267 


-6.890894 


-6.890894 




[£ = 0, n = 1] 


-71.604123 


-71.561383 


-71.604123 


-71.604093 


-71.604071 


[£ = 0,n = 2] 


-27.101595 


-26.875013 


-27.102023 


-27.101544 


-27.101509 


[£ = 1, n = 1] 


-50.549663 


-50.462831 


-50.549663 


-50.549664 


-50.549664 


[1 = 1, n = 2] 


-7.974057 


-6.549628 


-8.004096 


-7.974058 


-7.973692 


[1 = 2, n = 1] 


-28.555288 


-28.368573 


-28.555288 


-28.555288 


-28.555265 


[£ = 3,n = l] 


-6.890894 

6 = 1.5 fm 


-6.352976 


-6.890894 


-6.890894 


-6.890848 


[£ = 0, n = 1] 


-71.604133 


-71.241552 


-71.604133 


-71.604093 




= 0, n = 2] 


-27.101611 


-26.431344 


-27.393398 


-27.101544 




[£ = 1, n = 1] 


-50.549663 


-49.640614 


-50.549663 


-50.549664 




= 1, n = 2] 


-7.974057 


-7.578719 


-8.724761 


-7.974058 




[£ = 2, n = 1] 


-28.555288 


-27.396344 


-28.555288 


-28.555288 




[f = 3, n = 1] 


-6.890894 


-5.845401 


-6.890894 


-6.890894 





pendix |0 we list the actual equations we need to solve. 
In all our CCE calculations the physical vacuum \Q n t) is 
represented by the corresponding TLO n i wave function. 
The discrete spectrum of the hamiltonian is then deter- 
mined one state at a time, by solving the i?-set of equa- 
tions described above. In this approach one may hope to 
account exactly for the existence of the continuum part 
of the spectrum, by introducing the hole-functions \xni) 
and effectively mapping this continuum onto the bound- 
state part of the spectrum. 

We consider here two qualitatively different scenarios: 
The first scenario is identical in spirit with the approach 
used in the recent implementations of the CCE 0, 
for spin-isospin shell saturated nuclei: we adopt a one- 
basis realization of the CCE, where we use the remaining 
part of the TLO basis (after eliminating those state used 
to represent the physical vacuum, and possible all other 
lower energy states of the same symmetry) , to construct 
the desired hole- function \x-ni)- The functions in this set 
arc definitely orthogonal to each other, but they do not 
form a complete basis in the Hilbert space. However, the 
condition (4> n i\Xnij = must be satisfied, and indeed 
the above set of functions is complete in the subspace 
of all possible \xn,e)- However, a simple manipulation of 
the equations reveals that introducing the hole-function 
formalism in this picture does not amount to any new 



information. We are working within the confines of a TLO 
potential, and no structure outside this space can be built 
into the desired \^ n i) cigenfunction of the hamiltonian. 
We have effectively discretized the continuum, and the 
accuracy of the calculated eigenvalues and densities is 
generally close to the results in the Heisenberg picture. 

The ground state [£ = 0, n = 1] has a variational min- 
imum Eh for a value b = 1.4 of the oscillator parame- 
ter, and the CCE result in this representation is more 
sensitive to the choice of b than the previous result in 
the Heisenberg representation. This reflects the impor- 
tance of a good choice for the physical vacuum for the 
CCE: a true variational calculation is necessary to es- 
tablish the best possible choice of b. In reviewing the 
results, we note on the plus side that this version of the 
CCE fixes some of the erratic behavior of the density for 
the \l = 3, n — 1] state. Still, we notice that too much 
strength is pushed into the interior region, as the falloff 
of the tail is too sharp. On the other hand, the CCE 
misses completely the second state for I = 1. From ta- 
ble[]]we see that already at the physical vacuum level, the 
TLO representation is unable to provide an acceptable de- 
scription of the [£ = 1, n = 2] state: the vacuum energy 
Eh={i2} does not have a minimum for b G [1.3, 1.5], and 
the subsequent calculation of the hole- function as a linear 
superposition of TLO functions is unable to remedy this 
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FIG. 1: (Color online) Density probabilities p(r) of WS bound states, derived using the "exact" shooting method (line), 
Heisenberg representation in a TiO basis (boxes), CCE using a single TtO basis (empty circles), and CCE using a combination 
of a HO and C basis (filled circles). Where appropriate, the values of the employed parameters were b = 1.4 fm, and 7 = 3 fm~ . 
All discrete basis featured N g = 68 functions, corresponding to the identical number of collocation point for the Gauss-Hermite 
or Gauss-Laguerre quadrature method, respectively. 



deficiency. Once again, the wrong asymptotic behavior 
of the TLO functions results into a inadequate representa- 
tion of the continuum, and thus the serious shortcomings 
of the CCE in this representation. 

In order to establish that this is indeed the case, we 
consider a second (two-basis) scenario in the CCE. We 
still represent the physical vacuum as a simple 7iO n e 
wave function, but we expand the hole-functions using 
Laguerre (£) functions as our second complete basis. Un- 
like the HO functions which falloff like cxp[-r 2 /(26 2 )], 
the asymptotic behavior of the C functions has the form 
exp(— 7r), similarly to the bound-state coulomb func- 
tions. Properties of the C functions are summarized in 
Appendix [D] and for the purpose of this two-set rep- 
resentation of the CCE, we are using a Gauss-Laguerre 
quadrature scheme with N g = 68 collocation points. 

It is interesting to remark that a CCE calculation done 
entirely in L space (similar to the CCE calculation in the 
TCO space described above), results in good eigenvalues 
for values of 7 = 3^4, but fails to generate reason- 
able densities: if the TLO functions were falling off too 
quickly with r, the C functions are spread out over large 
distances and have little usable strength in the interior 
region which can be used to build the density. When 
the two basis are combined, however, the CCE results 
in the two-basis representation literally collapse onto the 
exact results, for both the bound state energies and den- 
sities. The numerical results show very little sensitivity 
with the value of the oscillator parameter b, for a fixed 
value of 7 = 3fm _1 : numerical differences for the ener- 



gies calculated for various values of b appear beyond the 
six significant digits presented here. The fact that we are 
working with a finite number (N g = 68) of basis states 
is evident as the CCE result is not the same as the exact 
results, but the error is less than 0.005% . 

The CCE calculation has no problem producing accu- 
rate energies and densities for the shallow bound states 
(E > — 10 MeV), and the CCE successfully corrects for 
the lack of a viable TCO representation in the case of the 
physical vacuum for the [£ = 1, n = 2] state. We would 
argue that a good physical vacuum representation is still 
important, but a good representation of the continuum is 
indeed crucial for an accurate calculation of the bound- 
states spectrum of the hamiltonian. 

We conclude this analysis by considering the normal- 
ization issues we have alluded to at the end of sec- 
tion IVII we remind the reader that the CCE wave func- 
tion is not normalized in the usual sense, but by requiring 
($|^) = 1. Therefore, the observables' calculation must 
implicitly "correct" for this artifact, and "restore" the 
wave function normalization, sec Eq. I|81|) . In table [HI 
we collect the results regarding the normalization of the 
solution and the normalization of the probability 
density p(r): for an exact solution of the Schrodingcr 
equation, the density should come out automatically nor- 
malized to unity, even though \^) is not. Conversely, the 
departure from unity in the p(r) normalization, is a test 
for the quality of the CCE solution. Therefore, it comes 
to no surprise that the largest discrepancies arise in the 
normalization of the [t = 3, n = 1] and, particularly, 
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TABLE II: Restoration of normalization of WS eigenfunctions 
in the CCE: using a pure TiO basis, and a combination of a 
TLO and C basis, respectively. All calculations correspond to 
the values of the parameters b = 1.4 fm, and 7 = 3 fen , 
respectively. 





state 




<*|*)-HO 


Pho 


(*|*)wO+£ 


PHO+C 


[1 


= 0, 


n — 


1] 


1.000979 


1.0005641 


1.000415 


1.000000 


[t 


= 0, 


n — 


2] 


1.000011 


0.9954375 


1.004602 


1.000000 


\t 


= 1, 


n — 


1] 


1.000157 


0.9988504 


1.001309 


1.000000 


V 


= 1, 


n — 


2] 


1.006852 


0.9323478 


1.081757 


1.000000 


[i 


= 2, 


n — 


1] 


1.001492 


0.9985834 


1.002913 


1.000000 




= 3, 


n — 


1] 


1.000184 


0.9823059 


1.018200 


1.000000 


[t 


= 1 


, n ■■ 




2] states, 


in the one- 


-basis TCO representa- 



tion. The two-basis representation of the CCE on the 
other hand, performs very well: the density normaliza- 
tion is correct to 11 significant figures. 



IX. DISCUSSION AND OUTLOOK 

In this paper we have discussed the fundamental as- 
pects of a new formulation of the coupled-cluster expan- 
sion in the continuum, and we have illustrated this ap- 
proach by considering the simplest (A = 1) dimensional 
realization of the many-body problem. For this simple 
scenario we have derived CCE equations which can be 
solved either in momentum or coordinate space. We have 
discussed the calculation of the eigenvalue spectrum of 
the hamiltonian, as well as the calculation of expectation 
values of arbitrary operators. 

The momentum space approach seems to be the most 
promising one for future efficient algorithms, with the 
drawback that an explicit knowledge of the continuum 
basis states appears to be mandatory: One can use 
Eqs. (|55|l and H56fl to design an iterative, albeit nonlinear, 
approach for the direct calculation of the energy, without 
an explicit reference to the S\ correlations. This may 
allow for a solution based on a nonlinear Monte Carlo 
approach, contingent upon the favorable resolution of se- 
rious concerns related to the numerical convergence and 
stability, and the feasibility of a practical numerical im- 
plementation. This discussion is tied in with the conver- 
gence and stability issues of the iterative method outlined 
in Section lYlII following Eq. (|88|l . Work is currently on- 
going to solve the technical issues of this numerical ap- 
proach. 

The coordinate space approach is also promising, and 
likely to offer immediate results. It offers the familiar 
perspective of a coordinate space implementation, with 
the added advantage of a lossless account of the contin- 
uum part of the single-energy spectrum. This is realized 
via a small set of /joZe-functions, labelled after the states 
which are occupied in the physical vacuum. However, 



a Monte Carlo approach to solving this set of equations 
seems difficult to envision at this time: The calculation 
of the ZioZe-functions employs an expansion in a known 
(and hopefully finite) basis of functions, with the expan- 
sion coefficients obtained as the numerical solution of a 
system of linear equations derived from Eq. (|7U|) . 

The stage is now set: We know that in the one-body 
case, we can derive consistent CCE equations. We will 
derive next the CCE equations for the deuteron (^4 = 2) 
problem. In the CCE formalism, one has an interme- 
diate level of approximation short of the exact solution, 
namely the Si and S2 levels, respectively. We will solve 
these equations using the realistic Argonne vig poten- 
tial [f|, compare the approximate and exact CCE so- 
lutions, and compare with results obtained using other 
exact methods. In the end, the deuteron problem is ex- 
tremely simple outside the CCE formalism, since it can 
be formulated as an one-body problem in terms of the 
relative coordinate of the system of two nucleons. 
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APPENDIX A: (A = 1) COMMUTATORS 

We find useful to write the following anti-commutation 
properties 

{ C pu a j? 2 } = $PlP2 > { C a,^>ln 2 } = 0, (Al) 

{c+ ,aj. } = , {clMnJ = 5 nim2 . (A2) 

We list now the results used in deriving the CCE equa- 
tions for our toy single-particle quantum mechanics prob- 
lem: 

We obtain the relevant expectation values 

(*|ct C/3 |$) = 5 nanp , (A3) 
($|b„ b a p - a 4 C/3 |$) = 5p a pJ n/3nb ■ (A4) 

• [H,Si] 

= C U+ b «! <Wl) + (- 4i 5 ^) C P 
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Note that 

[ctc^a^btjld.) 

Therefore, the relevant expectation values are 



&ppPi&n a n\ j 



(A5) 
(A6) 



(*![,] I*) 

($|b„ 6 a p - a [,]|*> 
' [[H,Si], Si] 

— — A A- - at ut _ A A- - Kt 

— < J n a n 2 u pppi cl p' 2 lJ n 1 u n„rii u p,3P2 "jfi u n 2 ■ 

We have 

(*l [[,],]!*> = 0, (A7) 
($|b„ 6a#o [[,],]|$} (A8) 

^n ci ri2^Pi3Pl^P2Pa^ n l n b ^ n a n l^P{3P2^PlPa^'>T>2'n>b ' 



APPENDIX B: (4 = 1): EQUATIONS WE NEED 
TO SOLVE 



In general the integral over the continuum particle 
states must be viewed as an integral over the continuum 
basis states, plus a sum over the set of discrete states 
in the spectrum, which remain unoccupied in the physi- 
cal vacuum |$). For simplicity, in the main body of the 
paper we assumed that all particle states are located in 
the continuum and this explicit sum was left out. In this 
appendix we revisit this issue and list the equivalent of 
Eqs. (E21 ED E3 El EZ£9, for the calculation of the 
energy and Iplh correlations. We also list the equivalent 
of Eqs. d3 E3 E3 HSU! required for the calculation of 
observables in momentum and coordinate space. 



1. Correlations and energies 

• Eigenstate of the hamiltonian H : 

*> = \4>h) + E 5 «^«> + / S p'^p) • ( B1 ) 

• Hole- function : 

\ Xh ) = |*> - |&> - £ SUK> ■ (B2) 



• P-set equations : 

E = {<t> h \K\cl> h ) + Y,{<t>h\K\<l>n) S nh 

d 3 p 



(27T) 



H|^)^, (B3) 



b n \R\<t> h )+ (Mn\^ ni ) s nif 



I 



bn]!!]^) S p - lh 



E Sa 



ph 



d 3 pi 
(2tt) 3 

p -|H|^)+ (0p|H|0 ni )5„ lh 



(B4) 



(H-£)| X = 

- (H-E) [\<j> h } + J2 \<f>m)S nih 



2. Observables 



P-set complement 



Snh — 



n h 



1 ^ SmhSmh 



d 3 Pl 



<5p7i — Sph 



(27T) 

1 ^ ^ S ni hS rL1 



3 ^p! /i *^pi ft; 



d 3 Pl 



3 ^pi/t^pi 



o = («^|0|*> 



(27T) 



6„|0|*)-(^|0|*)5„ fc 



d 3 p 



(B5) 



(2tt) 3 
• i?-set equations : 

£ = (^|H|0 h ) + <&|H| X fc) 

+ ^(0 /i |H|0„)^, (B6) 

= (4> n \n\(t> h ) + <^»|H|xfc> 

+ E (0„|H|0 ni )5 ril/l , (B7) 



(B8) 



(B9) 



(BIO) 



(Bll) 
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• i?-sct complement : 



where 



l-i<Xfc|H|*) 



(B12) 



- E s nih -(<p ni \m 



(Xh\ 



l-~<Xh|H|*> 



(B13) 



£ S nih -^ ni |H|*) 



6 = (0 h |O|*) + {x h |O|*)--g(xfc|H|*)(0 fc |O|*) 
+ £ [(4> n \0\V) - | (^„|H|*)(^|0|*)] . (B14) 



APPENDIX C: i?-SET PRACTICAL 
IMPLEMENTATION 

Consider two complete sets of functions, \<j) n ) and \tp n ). 
We model the physical vacuum using elements of the first 
set, while the second set is used to determine the corre- 
sponding hole-function. In general, we expect \4> n ) to 
offer a good description for the interior part of the wave 
function, while \tp n ) will provide a good description for 
the tail of the wave function. In practice, we use har- 
monic oscillator and Laguerre wave functions as the set 
\<p n ) and the set \ip n ), respectively. Note that a differ- 
ent viable candidate for the second basis may be the set 
of transformed harmonic oscillator wave functions intro- 
duced in |4lj . functions which are obtained via a local- 
scale point transformation of the spherical harmonic os- 
cillator basis. At this time we find the Laguerre functions 
simpler to handle. Properties of the Laguerre functions 
are reviewed in Appendix ID1 

We list now the necessary equations for the ground 
state. Extension of the formalism to the excited states 
calculation follows without difficulty. For the ground 
state, we take |<J>) = \<j>\)- We write the hole-function 
as an expansion in the second set of functions: 



\Xi) 



E C„ |V>n 



(CI) 



with |^„) = \ip„) - |</>i) (0i |?/> n ). Since \ip n ) are not or- 
thogonal to |</>i), the overlap substraction is necessary to 
fulfill the condition (xi \<j>x) = 0. The price we have to pay 
is that the set of functions \4> n ) are no longer orthogonal 
to each other. The expansion coefficients c n are found as 
the solution of the linear system of equations i|70|) : 



E 



(ip m \H\ip v 



E {4> m \4>n 



$m|H|0!> , 

(C2) 



(ViJHK) = (VUH|VV>)+£l<V'm|<M<4>l|V>n) (C3) 

- (V m |0i)(0i|H|V„) - (^ m |H|^i)<0i|V„> , 

and 

(V> m |-0n) = Smn ~ (ipm 1 01 ) (01 1 Ipn) ■ (C4) 

We calculate the ground state energy by using Eq. (|66[) : 
E = (0i|H|0i) + E c » (^i|H|-0„> , (C5) 

where 

(<f>i\K\iJ n ) = (^ilHlV^-^i^ilVn)- (C6) 

For the calculation of observables we also need the 
modified hole-function, |xi). Once again we perform an 
expansion in terms of the set \ip n )- We have 



\X 



1) = E dn i^") ' 



(C7) 



where the expansion coefficients are found by solving the 
linear system of equations (|79|) : 



E rf « (V>n|VO+ p (Xl|V>m) WV,|H|¥) = US 



E 



(C8) 



Subsequently, observables are calculated using Eq. (jHUJ). 
For local operators, it is sufficient to calculate the one- 
body density 

p(f) = (*|<y(f-fi|*)/(*|*) (C9) 

= ^(Wf)[l-i<Xi|H|*)] + *T(f)*(f) • 

Then, the expectation value of an arbitrary one-body 
local operator O can be calculated as 

<tf|0(r=)|#) = /d 3 n O(fi) (#|5(f-fi|*) . (CIO) 



APPENDIX D: SINGLE-PARTICLE WAVE 
FUNCTIONS 

We briefly review here the properties of the single- 
particle wave functions employed in the numerical part 
of this paper. We follow the conventions of Ref. [42^ . 

• the harmonic oscillator functions have the form 
HO nl {x) ~ x l+1 e- x2 ' 2 4 +1 'V) , (Dl) 
with x — r /b, and obey the orthogonality condition 



a <un < < ^ (2» + 2l + l)!!^ 

(D2) 
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Here {L^(x) ; n = 0, 1, . . .} arc the associated La- 
guerre polynomials. 

For the calculation of the kinetic energy matrix el- 
ement we will need the identity 



d 2 1(1+1) 
dr 2 

1 



b 2 



2(2n- 



3, , 



(D3) 



no ne (x) , 



and the nonzero elements of the (n£\x\m£) matrix 

3 



(n£\x 2 \n£) = 2n + £ + 



(n-1 £|z 2 K) 



n(n + £+ -) 



(D4) 
(D5) 



• the Laguerre functions have the form 

C^\x) ~ *"V-' a ^(a) , (D6) 
where x = jr, and obey the orthogonality condition 



' dx £^ a (x) £%+»{x) = {n + 2t + 2)l 6 n 

nl 



(D7) 



For the calculation of the kinetic energy matrix el- 
ement we need the identity 



d 2 _ 1(1+1) 
dx 2 



.r- 



1 n+£+ 1 n 
A x x 2 

n + 2£ + 2 „ 2 £+2 



^ +2 (*) 



(D8) 



(D9) 



The above equation involves the unnormalized La- 
guerre functions. If the normalized Laguerre func- 
tions are desired, then the coefficient of C?^-\ must 
be modified accordingly. 

The following identities pertain to the matrix ele- 
ment calculation of (n£\x~ p \m£) . with p= 1,2. We 
have 



dx x M+ " e~ x Ll l+2 {x)L^ +z {x) (D10) 



i{n,m} 



E 



{21 + 1 + c) ! 



and 



dx x™ e- l%+\z)I%»{x) 

min{n,m} 

(n + l-c)(m+l-c) 

c=0 



(Dll) 



(2£ + c) ! 
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